When does cyclic dominance lead to stable spiral waves? 
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Species diversity in ecosystems is often accompanied by characteristic spatio-temporal patterns. 
Here, we consider a generic two-dimensional population model and study the spiraling patterns 
arising from the combined effects of cyclic dominance of three species, mutation, pair-exchange and 
individual hopping. The dynamics is characterized by nonlinear mobility and a Hopf bifurcation 
around which the system's four-phase state diagram is inferred from a complex Ginzburg-Landau 
equation derived using a perturbative multiscale expansion. While the dynamics is generally char- 
acterized by spiraling patterns, we show that spiral waves are stable in only one of the four phases. 
Furthermore, we characterize a phase where nonlinearity leads to the annihilation of spirals and to 
the spatially uniform dominance of each species in turn. Away from the Hopf bifurcation, when the 
coexistence fixed point is unstable, the spiraling patterns are also affected by the nonlinear diffusion. 
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In nature, organisms live in areas much larger than the 
distances they typically travel and thus interact with a 
finite number of individuals in their neighborhood. Space 
and mobility are therefore crucial ingredients to under- 
stand how populations evolve, and how ecosystems self- 
organize. They are often responsible for pattern forma- 
tion, whose origin in ecosystems has been a subject of 
intense research for decades [TJ [2]. In his pioneering 
work, Turing showed that pattern-forming instabilities 
can be caused by diffusion [3J. While Turing patterns 
have been found in ecology and biology pQ, the require- 
ments of Turing's theory (e.g. separation of scales in 
diffusivities) appear to be too restrictive to explain the 
formation of patterns in many ecosystems, see e.g. [1]. 

Another important problem concerns the mechanisms 
promoting the maintenance of biodiversity [5]. In this 
context, cyclic dominance has recently been proposed 
as an intriguing motif facilitating the coexistence of di- 
verse species in ecosystems. Examples of cyclic compe- 
tition between three species can be found in coral reef 
invertebrates, Californian lizards, and communities of 
E. coli [3J [7J . In the experiments of Ref. [B] , the cyclic 
competition of three bacterial strains on two-dimensional 
plates was shown to yield intricate patterns sustaining 
species coexistence. Such competition is metaphorically 
described by rock-paper-scissors (RPS) games, where 
"rock crushes scissors, scissors cut paper, and paper 
wraps rock" [Sj. While non-spatial RPS-like models often 
evolve toward extinction of all but one species in finite 
time [5], their spatial counterparts are generally charac- 
terized by the long coexistence of species and by the for- 
mation of complex spatio-temporal patterns [T0lfT4] . Re- 
cently, two-dimensional versions of a model introduced by 
May and Leonard [TS] have received much attention [TlT - 
IT4] . When mobility is implemented by pair-exchange 
among neighbors, species coexistence is long-lived and 
populations form non- Turing spiraling patterns below a 
certain mobility threshold, whereas biodiversity is lost 
when that threshold is exceeded 



In this Letter, we characterize the intricate patterns 
emerging from the dynamics of a generic model of a cycli- 
cally competing three-species population, and study how 
these patterns affect the maintenance of biodiversity in 
two dimensions. The basic evolutionary processes con- 
sidered are the most general form of cyclic dominance 
obtained by combining and unifying the interactions of 
Refs. HUES [HUB]. Ins P ired b y the experiments of [7J, 
the model is formulated at the metapopulation level [TTJ , 
and is characterized by a Hopf bifurcation as well as by 
nonlinear mobility. While spiraling patterns have often 
been observed numerically in related models [111 - 114] . we 
here demonstrate that nonlinearity and mobility can dis- 
rupt the stability of the ensuing spiral waves. Our main 
result is the state diagram derived from a controlled mul- 
tiscale expansion around the Hopf bifurcation's onset and 
characterized by three phases in which spiral waves are 
unstable, and by one phase where spiraling patterns are 
stable. This implies the existence of a region of the pa- 
rameter space where spiral waves annihilate and each 
species dominates the system in turn. 

The generic model of cyclic competition is defined as a 
periodic square lattice of L 2 patches (L being the linear 
size) labeled by a vector I — (I1J2) [IB]- Each patch 
has a limited carrying capacity, accommodating at most 
TV individuals, and consists of a well-mixed population 
of species Si, S2, S3 and empty spaces 0. Within each 
patch I, the population composition evolves according to 

S, + S !+1 A S, + S, + S l+1 ^ 2S 4 (1) 
S + 0A2S, SASfcti, (2) 

where the species index i S {1,2,3} is ordered cycli- 
cally such that S3+1 = Si and Si_i = S3. The 
reactions ([!]) describe the cyclic competition between 
the species: Sj dominates over Si+i while being dom- 
inated by Si_i. Here, we consider a generic form of 
cyclic competition by separating the zero-sum process of 
dominance-replacement (rate £), as studied in Ref. |12) . 
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from the dominance-removal selection process (rate a) 
of Refs. [TTJ UM- With reactions @, we assume that 
births (rate 0) occur independently of the cyclic compe- 
tition provided that space (0) is available, and that each 
species can mutate into another (rate /i). 

As biological movement is often nonlinear and driven 
by local population density [TH] , we here divorce hopping 
(rate 5d) from pair-exchanges (rate Se) between nearest- 
neighbor patches I and V , according to 



[ s i\iWi> 
[Si] j [5 i± i] j, 



MM* 



(3) 



where I and I' lie in 4-neighborhood. The processes ^ 
lead to nonlinear mobility (see ^ below) and allow us to 
distinguish the movement in crowded regions, where pair- 
exchange dominates, from mobility in dilute systems, 
where hopping is more likely. The metapopulation model 
|l|-([3]) is well-suited to capture stochastic effects via size 
expansion in the carrying capacity and allows a natural 
connection with its deterministic description i4, 201 121] . 

It has to be noted that most previous works considered 
lattice models with N = 1 and nearest-neighbor reactions 
Q-Q, while here these interactions occur on-site. Apart 
from these differences, the processes that we consider are 
similar to those of [TT] when £ = fi = and 8r> = Se, 
while some aspects of the system's properties with £ ^ 0, 
fi 7^ and Sn ^ Se have been investigated in [12], [16] 
and [H], respectively. 

When N — > oo, the leading-order term in the size ex- 
pansion yields mean field rate equations for the contin- 
uous species densities Sj = Ns t /N. Here, Ns ( is the 
number of Si's in one patch. With s = (si, S2, S3), 



dsj 
dt 



Sj[/3(1 - r) - asi-t + ((s i+1 - s;_i)] 
+H(si-! + s i+1 - 2si) = Fi(s), (4) 



where r = si + S2 + S3 is the total density. Eqs. Q 
admit a coexistence fixed point s* — 3 ^ +cr (1, 1, 1)- In 
the presence of mutations, s* is an asymptotically stable 
focus when jj, > /i# = sTgf^) > while there is a supercrit- 
ical Hopf bifurcation (HB) [TO] at fj, — /i# and a stable 
limit cycle of frequency luh ~ ^2(3^^ when /1 < 
For later convenience, the departure from the HB point 
is measured by a parameter e defined by /1 = /.in — \t 2 ■ 
In stark contrast, when fj, = (no mutations), the coex- 
istence state s* is never asymptotically stable. Instead, 
solutions of Q are either heteroclinic cycles (/i = and 
er > 0) [TS] or nested neutrally stable periodic orbits 
(/i = a = 0) [8]. In either case, finite-size fluctuations 
cause the rapid extinction of two of the three species in 
a non-spatial setting [9]. 

When spatial dependence is taken into account in the 
limit L — > 00 and lattice spacing —¥ 0, the spatial coordi- 
nate x = l/L becomes continuous. The densities depend 




FIG. 1: Reactive steady states in stochastic simulations of 
reactions Here, L 2 = 128 2 , N = 64, /3 = a = 5 D = 

= 1, /Lt = 0.02 < fi H = 0.042 (e = 0.25) and, from left to 
right, £ = (1.8, 1.2,0.6,0). Each pixel describes a patch with 
normalized RGB representation (red, green, blue) = (si, S2, S3) 
of its state. The right-most panel shows an oscillatory ho- 
mogeneous state in which each of the species dominates the 
whole population in turn (see Fig. [3] for time evolution). Ini- 
tially s w s* with small random perturbations, see [18] . 



on space and time, Si = Si(x,i), and obey 
d tSl = Fi(s) + S D A Sl + (S D - S E ) {s l Ar - rA Sl ) , (5) 



where A = d£ + d^. 2 and the nonlinear diffusive 
terms (s^Ar — rAs^) arise from the divorce between 
pair-exchange and hopping. With our metapopulation 
approach, these partial differential equation (PDEs) can 
be derived in the continuum limit at the lowest order 
of a size expansion in N of the Markov chain associated 
with the processes [4] [20]. 

We aim to unravel the combined influence of nonlin- 
earity, mobility and noise on the system's dynamics and 
the formation and stability of coherent patterns. To 
gain some insight, we report some typical lattice sim- 
ulations obtained in the regime where there is a limit 
cycle (fi < /iff). As shown in Fig. [l] this regime is 
characterized by spiraling patterns found in four different 
phases, whereas we have found no patterns when fi > /in 
(see [15]). We have checked that the PDEs ^ faithfully 
reproduce the behaviors obtained with lattice simulations 
as shown in Fig. [l] Our analysis is therefore based on ([5| , 
whose properties are conveniently discussed by perform- 
ing the linear transformation s — s* — > (u,v,w), with 
u = — (r + S3) /a/6, v = (S2 — Si)/-\/2 and w — r/yo. In 
these variables, the linear part of Q can be written in 
the Jordan normal form d t (u + iv) = (e 2 + iuj^iu + ?w ) 
and d t w = —j3w. 

To make analytical progress, we perform a space and 
time perturbation expansion in the parameter e around 
the onset of the HB [22 . For this, we introduce multiple 
scale coordinates T = e 2 t and X = ex with Ax = d\ 



f 1 

, and expand the densities in powers of e. This yields 



u(x,t) = J2^U {n) (t,T,X) 



(6) 



n=l 



and, similarly, v = ELi e "^ (n) and w = EiUi e"^ (n) , 
where the functions ,V^ n \W^ are of order 0(1). 
Substituting Q into ^ and, using the definition of 
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FIG. 2: Upper: System's state diagram around the HB's on- 
set with contours of c = (cai,cei,cbs) f° r P = 1- We dis- 
tinguish four phases: spiral waves are unstable in AI, EI and 
SA, whereas they are stable in BS (see text). Lower: Typical 
snapshots from the PDE |5| in phases AI, EI, BS, SA from 
left to right (compare with Fig. [l] same parameters used). 

(u,v,w), we obtain a hierarchy of PDEs and analyze 
them at each order of e. To obtain a sensible expan- 
sion all secular terms are removed. Since the variables u 
and v are decoupled from w at linear order, one writes 
[/(i) + iv^) = A{T,X)e 2UJH \ where A is the complex 
modulation amplitude. The decoupled equations for w 
give W( 1] = and oc \A\ 2 , which is the lead- 

ing term in the equation for the center manifold |23j . 
The first secular terms arise at order C(e 3 ) and their 
removal yields the complex Ginzburg-Landau equation 
(CGLE) [53] with a real diffusion coefficient S 

d T A = SA x A + A-(l + ic)\A\ 2 A, (7) 

where S = 3l3S 3 E p+° SD and A has been rescaled by a con- 
stant to give 

12C(6/3^ ( T)(a + C) + ^(24/?-a) 

c = ^ . (8) 

3V3a(6/3 + a)(a + 2() 

We emphasize that the CGLE Q has here been derived 
perturbatively and describes the system's dynamics to or- 
der e near the HB's onset. This treatment, therefore, 
differs from that of Refs. Q3HH US], where a CGLE is 
obtained by treating heteroclinic cycles as limit cycles. 

According to the CGLE ([7]), the movement in the vicin- 
ity of the HB is described by linear diffusion, with an ef- 
fective diffusion constant S depending on So and Se 
When reproduction dominates over selection {(3 ^> a), 
the lack of empty spaces leads to prevalence of pair- 
exchanges (5 — > Se), while in the opposite case (/3 <C er), 




FIG. 3: Typical time evolution of the stochastic system in the 
SA phase. Same parameters and initial conditions as in Fig. [I] 
with C = 0. Upper: spiral annihilation at different stages, for 
time t — (234,310,386) from left to right. Lower: the oscil- 
latory dominance of each species at t — (955, 967, 980) after 
relaxation into the homogeneous state (no species extinction) . 



movement occurs mostly via hopping (5 — > <5d). As the 
effective linear diffusive term in ^ affects only the pat- 
tern's size, S can always be rescaled to 1 via x — > x/y/5. 

The system state diagram (near the HB) can be in- 
ferred from ([7]) and ^ by referring to the well-known 
properties of the two-dimensional CGLE [24] and is char- 
acterized by four phases with three critical values of c, 
as illustrated in Fig. [2j In the "spiral annihilation" (SA) 
phase, when < c < cbs, the dynamics is characterized 
by unstable spiraling patterns that collide and vanish. In 
the "bound state" (BS) regime cbs < c < cei, pairs of 
stable spirals are formed and coevolve, with their prop- 
erties described by the CGLE When cei < c < cai, 
the spirals become convectively unstable due to the Eck- 
haus instability (EI) which limits their size and distorts 
their shape. It is noteworthy that EI has been reported 
in [12] for a model without mutations (fi = 0). Finally, 
there is the "absolute instability" (AI) of spiral waves 
when cai < c, where there are no coherent patterns as 
the cores are not able to sustain spiral arms. With the ex- 
plicit values cbs ~ 0.845, cei ~ 1-25 and cai ~ 1-75 [2"4] . 
one obtains the system's state diagram in the a — C plane 
as shown in Fig. [2] This state diagram sheds light on 
the results of Fig. jTlwhere the values £ = (1.8, 1.2, 0.6, 0) 
correspond to c = (1.9, 1.5, 1.0, 0.6), which lie in the four 
phases AI, EI, BS and SA respectively. A description of 
the evolution in each phase can be found in the accom- 
panying movies [18]. The phase SA (see Fig. |3|, char- 
acterized by the annihilation of all spiraling patterns is 
particularly interesting since it is the only possible phase 
near the HB's onset when ( = 0, i.e. for the models 
of [TTJ [H] supplemented by mutations. In this phase, 
spiral annihilation leads to a spatially-homogeneous os- 
cillating state dominated in turn by each species, with- 
out any of them going extinct, as described by Q. This 
deterministic phenomenon is driven by nonlinearity (dif- 
ferent from the EI) and not by demographic noise. In 
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FIG. 4: Influence of mobility on spiraling patterns: 
typical snapshots in lattice simulations for (8d,Se) = 
(0.05,0.05), (0.20,0.05) from left to right respectively. Other 
parameters are: L 2 = 128 2 , TV = 1024, = cr = 1, C = 0.1 
and (i — 10 -6 <C /in = 0.042. Geometrically ordered initial 
conditions, see movies [18] for full details. 



bifurcation, are robust and independent of the mobil- 
ity rates. This is in stark contrast with the results of 
Refs. [HI [T2J, [14], where spiraling patterns and spatial 
uniformity were respectively found at low and high mo- 
bility, and may explain why spiraling patterns turn out to 
be elusive in the microbial experiments of Refs. [SJ[7] • Ad- 
ditionally, regardless of internal noise, we show that non- 
linear diffusion causes far-field break-up of spiral waves 
away from the Hopf bifurcation's onset. 

The authors acknowledge discussions with Tobias 
Galla at an early stage of this project. BS is grateful 
for support of an EPSRC studentship. 



the regime c <C cbs, it typically occurs on a short time 
scale, as illustrated in Fig. [3] This is markedly different 
from the loss of spiraling patterns driven by noise after a 
time growing exponentially with the system size as found 

in [nun]- 

Our analysis in terms of the CGLE ^ relies on a per- 
turbative treatment around the onset of HB where e«l, 
but is found to faithfully describe the system's properties 
far from the HB. For instance, when /3 = a = 1 and £ = 
(/j,h = 0.042), the system is still in the SA phase even 
for [i = 0.02 (e = 0.25) as predicted by our theory (see 
Figs. [T]and[3]). We have also checked that our analysis 
is robust against simultaneous random perturbations (up 
to ±5%) of all the reaction rates Q-Q [20]. As shown 
in Figs, [I] and [2j the PDEs ^ describe perfectly the 
stochastic metapopulation model when N ^> 1 and, in 
practice, are still accurate when N > 16 for any nonzero 
mobility. Furthermore, when N = 2 and the mobility 
rates are sufficiently high, as discussed in [II] . the state 
diagram of Fig. [2] is still valid [15] , 

Away from the HB's onset (e.g. for /i = 10~ 6 ) and 
when 8d ^ 5e, the stability of spirals is affected by the 
nonlinear diffusive terms of ([5]). In fact, when noise is 
negligible (N 3> 1), one observes a far-field break-up of 
the spiraling patterns caused by nonlinear mobility, as 
shown in Fig. [1] and [18] . 

In summary, we have investigated the stability of spi- 
raling patterns in a generic three-species model whose 
evolution results from the combined effects of cyclic dom- 
inance, mutation and biologically-motivated nonlinear 
mobility. Inspired by recent experiments [7] , we have de- 
veloped a metapopulation description and analyzed the 
dynamics in terms of a CGLE derived via a multiscale 
expansion around the Hopf bifurcation's onset. We have 
thus obtained the system's state diagram, which is char- 
acterized by four phases, with only one capable of sup- 
porting stable spiraling patterns. The instabilities in the 
three other phases are not driven by noise. In particular, 
we have identified a phase (SA) where spirals annihilate, 
leading to spatially uniform dominance of each species 
in turn. Importantly, these behaviors, which arise in 
a wide region of the parameter space around the Hopf 
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